{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://hilpisch.com/tpq_logo.png\" alt=\"The Python Quants\" width=\"35%\" align=\"right\" border=\"0\"><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Listed Volatility and Variance Derivatives\n",
    "\n",
    "**Wiley Finance (2017)**\n",
    "\n",
    "Dr. Yves J. Hilpisch | The Python Quants GmbH\n",
    "\n",
    "http://tpq.io | [@dyjh](http://twitter.com/dyjh) | http://books.tpq.io\n",
    "\n",
    "<img src=\"http://hilpisch.com/../images/lvvd_cover.png\" alt=\"Derivatives Analytics with Python\" width=\"30%\" align=\"left\" border=\"0\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Valuing Volatility Derivatives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chapter illustrates the valuation of volatility futures and options according to Gruenbichler and Longstaff (1996), also abbreviated in the following by GL96. They derive a semi-analytical (\"closed\") pricing formula for European volatility call options which is as easy to use as the famous Black-Scholes-Merton formula for equity options pricing. They model volatility directly and make the assumption that volatility follows a square-root diffusion process. The model is quite simple and parsimonious such that it lends itself pretty well to serve as a starting point.\n",
    "\n",
    "This chapter introduces the financial model, the futures and option pricing formulas of Gruenbichler and Longstaff as well as a discretization of the model for Monte Carlo simulation purposes. For both, the formulas and the Monte Carlo simulation approach, Python implementations are presented. Finally, the chapter shows in detail how to calibrate the GL96 model to market quotes for European call options on the VSTOXX volatility index."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Valuation Framework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gruenbichler and Longstaff (1996) model the volatility process (e.g. the process of a volatility index) in direct fashion by a square-root diffusion or CIR process &mdash; named after the authors John C. Cox, Jonathan E. Ingersoll, and Stephen A. Ross, who first introduced this type of stochastic process to finance; cf. Cox et al. (1985). The stochastic differential equation describing the evolution of the volatility over time takes on the form\n",
    "\n",
    "$$\n",
    "dv_t = \\kappa (\\theta - v_t)dt + \\sigma \\sqrt{v} dZ_t\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where\n",
    "\n",
    "* $v_t$ is the time $t$ value of the volatility index, for example the VSTOXX\n",
    "* $\\theta$ is the long-run mean of the volatility index\n",
    "* $\\kappa$ is the rate at which $v_t$ reverts to $\\theta$\n",
    "* $\\sigma$ is the volatility of the volatility (\"vol-vol\")\n",
    "* $\\theta, \\kappa, \\sigma$ are assumed to be constant and positive\n",
    "* $Z_t$ is a standard Brownian motion\n",
    "\n",
    "This process is known to exhibit convenient and realistic features for volatility modeling, like positivity and mean-reversion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the numerical examples in this chapter we use a parametrization as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "import math\n",
    "## model parameters\n",
    "v0 = 17.5  # initial level of volatility index\n",
    "# kappa = 1  # high speed of mean reversion\n",
    "kappa = 0.1  # slow speed of mean reversion\n",
    "theta = 20.0  # long-term index level\n",
    "sigma = 2.0  # volatility of volatility\n",
    "zeta = 0.0  # factor of the expected volatility risk premium\n",
    "r = 0.01  # risk-free short rate\n",
    "\n",
    "## option parameters\n",
    "K = 20.0  # strike\n",
    "T = 1.0  # time horizon in year fractions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Futures Pricing Formula"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Denote by $F(v,T)$ the futures price of a futures contract on the volatility index $v$ with maturity $T$. The pricing of such a contract is done by taking expectations of the index's value at maturity $F(v,T)=E(v_T)$. One obtains the following formula for the futures price:\n",
    "\n",
    "$$\n",
    "F(v_0,T)= \\left( 1-e^{-\\kappa T}  \\right) \\cdot \\theta + e^{-\\kappa T} \\cdot v_0\n",
    "$$\n",
    "\n",
    "It is an exponentially weighted average of the long-run mean and the current value of the volatility index. The function ``futures_price`` implements this formula in Python, now for a general volatility risk premium $\\zeta \\geq 0$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sys.path.append('./scripts')\n",
    "import srd_functions as srd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srd.futures_price??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Application of the formula and function are straightforward. However, we have to first import the function from the module in which it is stored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "srd.futures_price(v0, kappa, theta, zeta, T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We obtain a futures term structure by calculating the futures prices for different maturities (see the following figure)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from pylab import mpl, plt\n",
    "plt.style.use('seaborn')\n",
    "mpl.rcParams['font.family'] = 'serif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "maturities = np.linspace(0, 2, 24)\n",
    "futures_prices = [srd.futures_price(v0, kappa, theta, zeta, T)\n",
    "                 for T in maturities]\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(maturities, futures_prices)\n",
    "plt.xlabel('maturity')\n",
    "plt.ylabel('futures price');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Futures prices from today to a maturity of 24 months."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Option Pricing Formula"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A European call option on the underlying $V$ has a payoff function $\\max(V_T-K,0)$ at maturity $T$ where $K$ is the strike price of the option. Gruenbichler and Longstaff (1996) derive the following pricing formula for such a call option"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{eqnarray*}\n",
    "    C(v_0,K,T) &=& e^{-rT} \\cdot e^{- \\beta T} \\cdot V_0 \\cdot \\mathbf{Q} \\left(\\gamma \\cdot K; \\nu + 4, \\lambda \\right)\\\\\n",
    "    &+& e^{-rT} \\cdot \\left(\\frac{\\alpha}{\\beta}\\right) \\cdot \\left(1-e^{- \\beta T}\\right) \\cdot \\mathbf{Q}\\left(\\gamma \\cdot K; \\nu+2,  \\lambda \\right) \\\\\n",
    "    &-& e^{-rT} \\cdot K \\cdot \\mathbf{Q}\\left(\\gamma \\cdot K; \\nu, \\lambda \\right)\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "    \\alpha &=& \\kappa  \\theta \\\\\n",
    "    \\beta &=& \\kappa + \\zeta \\\\\n",
    "    \\gamma &=& \\frac{4 \\beta}{\\sigma^2 \\left(1-e^{- \\beta T}\\right)} \\\\\n",
    "    \\nu &=& \\frac{4 \\alpha}{\\sigma^2} \\\\\n",
    "    \\lambda &=& \\gamma \\cdot e^{- \\beta T} \\cdot v_0\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and $e^{-rT}$ as the discount factor for a fixed short rate $r$. The parameter $\\zeta$ denotes as before the expected premium for volatility risk. In the following, we assume $\\zeta = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{Q}(\\cdot; \\nu, \\lambda)$ is the complementary non-central $\\chi^2$ distribution with $\\nu$ degrees of freedom and non-centrality parameter $\\lambda$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our Python implementation ``srd_functions.py`` (see the appendix for the complete module), we use the function ``ncx2.cdf()`` from the ``scipy.stats`` sub-library for the non-central $\\chi^2$ distribution. We implement two different functions, a helper function and the pricing function itself:\n",
    "\n",
    "* ``cx()``: this function returns values for the complementary distribution of the non-central chi-squared density\n",
    "* ``call_price()``: the implementation of the valuation formula for European calls on the volatility index\n",
    "\n",
    "The function ``cx()`` is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srd.cx??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function ``cx()`` is used in the valuation function ``call_price()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srd.call_price??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, being equipped with such a Python function makes applying the pricing formula straightforward. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "srd.call_price(v0, kappa, theta, sigma, zeta, T, r, K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us calculate European call option prices over a wider range of strikes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "option_values = []\n",
    "strikes = np.linspace(15, 25)\n",
    "option_values = [srd.call_price(v0, kappa, theta, sigma, zeta, T, r, k)\n",
    "                 for k in strikes]\n",
    "data = pd.DataFrame(option_values, index=strikes,\n",
    "                    columns=['call values', ])\n",
    "data.index.name = 'strike'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the following figure shows, the European call option prices behave comparable to equity options: the higher the strike price the lower the option present value ceteris paribus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data.plot(figsize=(10, 6));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">European call option prices for different volatility strikes in GL96 model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we solve the pricing problem for a European call option via Monte Carlo simulation, i.e. based on a large number of simulated volatility index paths. To this end, we have to discretize the stochastic differential equation to obtain a difference equation which we can evaluate numerically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the same numerical parameters as before, but we have to define some additional parameters which will determine the discretization interval and the number of simulated paths:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## simulation parameters\n",
    "M = 50  # time steps\n",
    "I = 20000  # number of MCS paths"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a large literature avaible about the discretization of the square-root diffusion (cf. Hilpisch (2015), ch. 10). In the following, we will implement an exact discretization scheme."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we divide the time interval $[0, T]$ into $M$ time intervals with $M+1$ discrete points in time: $t \\in \\{0, \\Delta t, 2 \\Delta t , \\ldots, T \\}$. We have:\n",
    "\n",
    "$$\n",
    "v_t = \\frac{\\sigma^2 (1- e^{- \\kappa \\Delta t})}{4 \\kappa} \\tilde{\\chi}^2_{\\nu,\\lambda} \\frac{4^{- \\kappa \\Delta t}}{\\sigma^2 (1- e^{- \\kappa \\Delta t})} \\cdot v_s\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $s = t - \\Delta t$.\n",
    "\n",
    "$\\tilde{\\chi}^2_{\\nu,\\lambda}$ denots a non-central chi-squared distributed random variable with\n",
    "\n",
    "$$\\nu = \\frac{4 \\theta \\kappa}{\\sigma^2}$$\n",
    "\n",
    "degrees of freedom and non centrality parameter\n",
    "\n",
    "$$\\lambda = \\frac{4 \\kappa e^{- \\kappa \\Delta t}}{\\sigma^2 (1- e^{- \\kappa \\Delta t})} \\cdot v_s$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it may be more convenient to sample a chi-squared random variable instead of non-central chi-squared one, we use the following equations:\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\nu>1 &:& \\tilde{\\chi}^2_{\\nu,\\lambda} = (z+ \\sqrt{\\lambda})^2 + {\\chi}^2_{\\nu-1} \\\\\n",
    "\\nu \\le 1 &:& \\tilde{\\chi}^2_{\\nu,\\lambda} = {\\chi}^2_{\\nu+2P}\n",
    "\\end{eqnarray*}\n",
    "\n",
    "where $z$ is an independent standard normally distributed random variable and $P$ is a Poisson distributed random variable with intensity $\\frac{\\lambda}{2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We again divide the implementation in Python in two functions:\n",
    "\n",
    "* ``generate_paths()`` generates simulated volatility level paths\n",
    "* ``call_estimator()`` simulates the volatility process and calculates the Monte Carlo estimator for the European call option value\n",
    "\n",
    "The function to generate simulated volatility paths takes on a form as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srd.generate_paths??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An application of this function might look as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "paths = srd.generate_paths(v0, kappa, theta, sigma, T, M, I)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "The following figure shows the first 10 simulated volatility paths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "pd.DataFrame(paths[:, :10]).plot(legend=False, figsize=(10, 6));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">10 simulated volatility paths for GL96 model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Monte Carlo estimator for a European call option is defined by\n",
    "\n",
    "$$\n",
    "C(v_0,K,T) = e^{-rT}\\frac{1}{I} \\sum_{i=1}^{I} \\max(v_T^i-K,0)\n",
    "$$\n",
    "\n",
    "with $\\max(v_T-K,0)$ being the payoff function of the European call option. Here, $v_T^i$ is the $i$-th simulated value for the volatility index at maturity. The function to calculate the Monte Carlo estimator for the value of the European call option is presented below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srd.call_estimator??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srd.call_estimator(v0, kappa, theta, sigma, T, r, K, M, I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, let us calculate option prices over a wider range of strikes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "estimates = []\n",
    "for k in strikes:\n",
    "    estimates.append(srd.call_estimator(v0, kappa, theta, sigma,\n",
    "                                    T, r, k, M, I))\n",
    "data['estimates'] = estimates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure compares the Monte Carlo estimator values given the above parametrization with the European call option values obtained from the formula by Gruenbichler and Longstaff (1996)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data.plot(style=['b', 'r.'], figsize=(10, 6));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Monte Carlo estimates vs. European call values from formula in GL96 model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automated Monte Carlo Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section illustrates how to automate Monte Carlo-based valuation tests. Even if semi-analytical option pricing formulae are available &mdash; like in the Gruenbichler and Longstaff (1996) model &mdash; one needs in general numerical methods as well. One of the most flexible and powerful methods when it comes to the valuation of more complex, exotic instruments is Monte Carlo simulation (MCS)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When implementing MCS algorithms it is generally advised to benchmark the results from these algorithms against results from other valuation methods which are known to deliver exact values. Therefore, this section implements MCS for the Gruenbichler and Longstaff (1996) model and benchmarks the results against the semi-analytical formula for European volatility call options as presented earlier in this chapter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous section, MCS is already introduced and some benchmark results are presented. However, this section enhances the analysis and adds functionalities for the structured storage and systematic analysis of valuation results. In this regard, the Python library pandas is used in combination with the storage capabilities of PyTables, the Python wrapper for the HDF5 database file format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Automated Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following test procedures are based on the Python script ``srd_simulation_analysis.py`` (see the appendix for the complete code). To automate the valuation for a significant number of options, a set of parameter values and lists has to be specified. First, some imports and the parametrization of the financial model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -n 10,23p scripts/srd_simulation_analysis.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameter values chosen are not too unrealistic, for example, for the VSTOXX volatility index. Second, the following simulation parameters are assumed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -n 25,38p scripts/srd_simulation_analysis.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The single parameters and list objects have the following meaning:\n",
    "\n",
    "* ``write`` is just a flag to indicate that the results should be stored on disk\n",
    "* ``var_red`` contains the configuration for the variance reduction techniques, i.e. moment matching of pseudo-random numbers and antithetic variates\n",
    "* ``steps_list`` is the list with the different numbers for the time steps to be used for the discretization\n",
    "* ``paths_list`` contains the different numbers for the MCS paths for the volatility process\n",
    "* ``runs`` is a parameter which says how often the same option is valued for a single parameter combination\n",
    "* ``SEED`` is a seed value for NumPy's pseudo-random number generator; it is fixed to make it easier to compare results\n",
    "* ``PY1`` is the performance yardstick for absolute differences\n",
    "* ``PY2`` is the performance yardstick for relative differences\n",
    "* ``maturity_list`` is the list with times-to-maturity\n",
    "* ``strike_list`` finally is the list with the different strike prices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For details with regard to such Monte Carlo simulation analyses and the respective MCS background information refer to the book Hilpisch (2018, ch. 12). For a detailed exposition of Monte Carlo methods implemented in Python see also Hilpisch (2015, ch. 10)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The script ``srd_simulation_analysis.py`` used in this section relies on the MCS function from the previous section. In addition, it uses the following function which implements two variance reduction techniques, moment matching and antithetic variates. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from srd_simulation_analysis import randoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "randoms??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main valuation and testing loop is implemented through the following code. Important passages are commented inline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -n 116,191p scripts/srd_simulation_analysis.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Storage Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The script ``srd_simulation_analysis.py`` from the previous sub-section generates a number of results and a lot of printed output during execution. However, one generally wants to store such results in a structured manner and analyze it maybe later on.  To this end, we use the pandas library which provides convenient data storage and analysis features. In combination with the PyTables library, simulation results are easily stored on disk and read from disk later on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider now the Python module ``srd_simulation_results.py`` with several functions to store and analyze the data generated through the other script (see the appendix for the complete script). We will go through this script part by part. Of course, there are some necesary imports and a filename for the database is also specified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -n 11,17p scripts/srd_simulation_results.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having a ``DataFrame`` object instantiated (with results already stored or not), the following function allows to add a set of valuation results to the ``DataFrame`` as an additional row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import srd_simulation_results as ssr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ssr.write_results??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once all the single simulation results have been added to the ``DataFrame`` ``sim_results``, this should be stored on disk in a HDFStore file object. This is what the following function accomplishes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ssr.write_to_database??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now turn to the analysis of valuation results. To this end, we work with the Python module ``srd_simulation_results.py`` (see the appendix for the code). The function ``print_results()`` prints simulation results stored in the HDFStore file object with name ``filename``. If not otherwise parametrized, the function prints the first 50 results stored in the database (if that many exist)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ssr.print_results??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us inspect selected valuation results from running the script ``srd_simulation_analysis.py``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ssr.print_results(idl=2, idh=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Time in Minutes\" and \"Time per Option\" values coincide since 60 options are valued per run.\n",
    "\n",
    "In principle, the results from these different parametrizations already illustrate that valuation accuracy is, in general, higher the higher the computational effort is (in this case the number of paths). The number of errors drops from 9 for 75,000 paths to 2 for 125,000 paths."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the second major function contained in the Python module is ``plot_error_ratio()`` and plots the error ratios against the computational effort, i.e. (number of simulation paths) x (number of time steps)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ssr.plot_error_ratio??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us generate an overview plot with the above function. The graphical output is shown in the following figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "filename = './data/simulation_results.h5'\n",
    "from srd_simulation_results import plot_error_ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plot_error_ratio(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Error ratios of Monte Carlo valuation vs. computational effort."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear regression illustrates graphically the statement that the higher the computational effort the higher in general the accuracy of the MCS valuations. The high variability of the error ratios stems from the fact that not only the computational effort influences valuation accuracy. For example, the use of variance reduction techniques also has a significant influence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This being said, one might be interested in analyzing the impact of the variance reduction techniques on the valuation accuracy. To this end, pandas provides powerful grouping methods for the ``DataFrame`` object. First, we open the database in which the simulation results are stored and read them into an in memory object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "h5 = pd.HDFStore(filename, 'r')\n",
    "data = h5['sim_results']\n",
    "h5.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation results are stored in 96 rows of the pandas ``DataFrame`` object. The method ``groupby`` allows analyses of such data sets in \"all directions and dimensions\". For example, one gets the average error ratio depending on the use of moment matching with the following code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data.groupby(['mo_match']).mean()[['error_ratio']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the simulations, moment matching obviously reduces the error ratio by a bit more than 2 percentage points from above 19% to less than 17%.\n",
    "\n",
    "The same analysis can be done with more data columns, both for the grouping and the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data.groupby(['mo_match', 'anti_paths']).mean()[['error_ratio', 'errors']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the results indicate that moment matching alone should be used. Adding antithetic paths increases the error ratio ceteris paribus.\n",
    "\n",
    "In similar fashion, the data can be grouped by time steps and paths &mdash; in analogy to the analysis above. The following figure presents the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df = data.groupby(['time_steps', 'paths']).mean()[['error_ratio', 'time_opt']]\n",
    "df.plot(x='time_opt', y='error_ratio', style='bo')\n",
    "plt.xlabel('time per option valuation')\n",
    "plt.ylabel('error ratio');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Error ratios of Monte Carlo valuation vs. time per option."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the general rule that accuracy increases with increasing computational effort is again supported by this specific analysis. Of course, other analyses can be done with the data generated and stored on disk."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Calibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Previous sections implement a *theoretical* valuation approach in that a parametric model is given and parameter values are simply *assumed*. However, any realistic valuation of volatility options has to incorporate at least some information from the market itself. This is generally done by taking the quotes of liquidly traded vanilla options as input for a so-called *model calibration*. During the calibration of a financial model, parameters are determined for that model which best replicate the observed option quotes. These parameters are then in turn used to value other, maybe more exotic, derivatives by the means of numerical methods like Monte Carlo simulation. Such a procedure then often is called a *market-based valuation approach*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section illustrates how to calibrate the financial model of Gruenbichler and Longstaff (1996) to option quotes for European volatility calls traded at Eurex. In the following, the semi-analytical formula as introduced previously is used to accomplish the model calibration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Option Quotes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The remainder of this section will work with option quotes of VSTOXX call options traded at Eurex which are all of European type. The data is as of 31. March 2014 and contains quotes for the options with maturities for each month from April 2014 to November 2014. Let us have a look at the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "path = './data/'\n",
    "h5 = pd.HDFStore(path + 'vstoxx_march_2014.h5', 'r')\n",
    "option_quotes = h5['vstoxx_options']\n",
    "option_quotes.info()\n",
    "n = 46244 - 25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "option_quotes = option_quotes.iloc[n:n + 98]\n",
    "option_quotes.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the data set comprises a total of 98 option quotes. Note that option quotes for Eurex volatility options have tick sizes of 5 Cent, i.e. 0.05 EUR. A detailed overview of the terms of these options can be found on the Web site http://www.eurexchange.com/exchange-de/produkte/vol/. In what follows, we assume that each option is written on a single unit of the underlying, i.e. 1 point of the VSTOXX index translates into 1 EUR."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Calibration Procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In simple terms, the problem of calibration is to find parameters for an option model such that observed market quotes of liquidly traded plain vanilla options are replicated as good as possible. To this end, one defines an error function that is to be minimized. Such a function could be the *Mean Squared Absolute Error (MSAE)*. The task is then to solve the problem\n",
    "\n",
    "\n",
    "$$\n",
    "\\min_{\\kappa, \\theta, \\sigma} \\frac{1}{N}\\sum_{n=1}^{N}\\left( C_{n}^{*} - C_{n}^{GL96}(\\kappa, \\theta, \\sigma)\\right)^{2}\n",
    "$$\n",
    "\n",
    "with the $C_{n}^{*}$ being the market or input prices and the $C_{n}^{GL96}$ being the model or output prices for the options $n=1,...,N$.\n",
    "\n",
    "For some constellations, it might be more appropriate to minimize the *Mean Squared Relative Error (MSRE)*:\n",
    "\n",
    "$$\n",
    "\\min_{\\kappa, \\theta, \\sigma} \\frac{1}{N}\\sum_{n=1}^{N}\\left( \\frac{C_{n}^{*} - C_{n}^{GL96}(\\kappa, \\theta, \\sigma)}{C_{n}^{*}}\\right)^{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Python script (the appendix for the complete code) that implements a procedure to calibrate the Gruenbichler-Longstaff option pricing model is explained in the following step-by-step.\n",
    "\n",
    "First, some imports and parameter specifications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -n 11,22p scripts/srd_model_calibration.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function ``read_select_quotes()`` reads and selects option quotes for the calibration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import srd_model_calibration as smc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smc.read_select_quotes??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the core functions of the calibration script are the following. On the one hand, ``valuation_function()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smc.valuation_function??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand, ``error_function()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smc.error_function??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function ``valuation_function()`` implements a valuation for the array of strikes under consideration and the provided parameter values. Function ``error_function()`` calculates either the MSAE or the MSRE &mdash; given the array of strikes and the parmeter vector ``p0``. Some penalties are also added depending on the single parameter values or their constellation. For details and a number of similar examples, refer to the book Hilpisch (2015, ch. 11).\n",
    "\n",
    "Then standard functions from the ``scipy.optimize`` library (cf. Hilpisch (2018), ch. 11) are used to implement both a *global* (i.e. via the ``brute()`` function) and a *local optimization* (i.e. via the ``fmin()`` function). This all is wrapped in the function ``model_calibration()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smc.model_calibration??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the results are visualized by plotting the model values against the market values and providing the absolute differences separately. This is accomplished by the function ``plot_calibration_results()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smc.plot_calibration_results??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Calibration Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start with a calibration run for the option quotes with the shorter maturity and using the *MSAE* as yardstick."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from srd_model_calibration import *\n",
    "option_data = read_select_quotes()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "option_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%time opt = model_calibration(option_data, rel=False, mat='2014-05-16')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the calibration results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plot_calibration_results(opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Calibration of option pricing model to option quotes for May 2014 maturity with MSAE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the calibration to the option quotes with a longer time-to-maturity and using the *MSRE* as yardstick."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%time opt = model_calibration(option_data, rel=True, mat='2014-07-18')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results of this calibration run are shown in the following figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plot_calibration_results(opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Calibration of option pricing model to option quotes for July 2014 maturity  with MSRE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Neither of the two fits is satisfactory for such small sets of option prices. The differences are often higher than the tick size of 5 Cent. Therefore, this model obviously does not paint a too realistic picture of the real world. This may be due to the underlying, i.e. the VSTOXX volatilty index, exhibiting more and/or other properties than the ones captured in the model of Gruenbichler and Longstaff (1996). For example, a well documented stylized fact is that time series of volatility indexes, like the VSTOXX, exhibit (positive) *jumps* of significant size and with significant probability. Something that is, for example, not captured in the pricing model used in this chapter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Volatility derivatives, like futures and options on a volatility index, and their arbitrage-free valuation are the focus of this chapter. The model of Gruenbichler and Longstaff (1996) is introduced since it serves pretty well as a benchmark case &mdash; due to its simplicity &mdash; for volatility modeling and volatility derivatives valuation. The model uses a square-root diffusion or CIR process to describe the evolution of volatility (indexes) over time. Futures and European call option pricing formulas are presented. Based on an exact discretization scheme, Monte Carlo simulation is described as an alternative numerical pricing method. In addition, Python code to implement automated Monte Carlo simulation studies is provided and used for such a study. Finally, the GL96 model is calibrated to market quotes for European call options on the VSTOXX volatility index. The results illustrate that the model might not be rich enough to account for market realities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Scripts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ``srd_functions.py``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat scripts/srd_functions.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ``srd_simulation_analysis.py``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat scripts/srd_simulation_analysis.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ``srd_simulation_results.py``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat scripts/srd_simulation_results.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ``srd_model_calibration.py``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat scripts/srd_model_calibration.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://hilpisch.com/tpq_logo.png\" alt=\"The Python Quants\" width=\"35%\" align=\"right\" border=\"0\"><br>\n",
    "\n",
    "<a href=\"http://tpq.io\" target=\"_blank\">http://tpq.io</a> | <a href=\"http://twitter.com/dyjh\" target=\"_blank\">@dyjh</a> | <a href=\"mailto:team@tpq.io\">team@tpq.io</a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
